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Abstract. In this work we propose a Bayesian framework for fully au- 
tomated image fusion and their joint segmentation. More specifically, 
we consider the case where we have observed images of the same object 
through different image processes or through different spectral bands. 
The objective of this work is then to propose a coherent approach to 
combine these data sets and obtain a segmented image which can be 
considered as the fusion result of these observations. 
The proposed approach is based on a Hidden Markov Modeling (HMM) 
of the images with common segmentation, or equivalently, with com- 
mon hidden classification label variables which are modeled by the Potts 
Markov Random Field. We propose an appropriate Markov Chain Monte 
Carlo (MCMC) algorithm to implement the method and show some sim- 
ulation results and applications. 
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1 Introduction 

Data fusion and multi-source information has become a very active area of re- 
search in many domains : industrial nondestructive testing and evaluation (PQ), 
industrial inspection (0), and medical imaging 3 4 5 1617] . For example in mag- 
netic resonance (MR) image segmentation, one can use multispectral images to 
accentuate the differencies in physical characteristics of the anatomical tissues. 
In this case we must fuse information of these multispectral images to obtain 
segmentation. 

The main problem is how to combine the information contents of different sets 
of data gi(r). Very often the data sets gi, and corresponding images /», do not 
represent the same quantities. A general model for these problems can be the 
following : 

g i (r) = [H i f i ](r)+e i (r), i = l,...,M (1) 
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where Hi are the functional operators of the measuring systems, or a registration 
operators if the observations have to be registered. We may note that estimating 
fi given each set of data gi is an inverse problem by itself. In this work we 
propose to reconstruct images fi and to construct a common segmentation at 
the same time. This fused segmentation will be presented by an image z. In this 
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Fig. 1. Examples of images for data fusion and joint segmentation, a) PD (Pro- 
ton Dennsity),Tl and T2-weighted slices of brain MR images, b) Tl-weighted, 
T2-weighted and Tl-weighted with contrast agent transversal slices of a 3D brain 
MR images with additive gaussian noise. 



paper we consider the case where the measuring data systems can be assumed 
almost perfect and the observations are registered, which means that we can 
write : 

9i(r)=f i (r)+e i (r), i = l,...,M (2) 

for r 6 R 2 . If we have, for example, multispectral images, then c/j's represent 
observations taken at different spectral bands. In a multimodal case, gi can 
represent CT and PET images. 

Our aim is then to obtain a common segmentation of N observations and to 
reconstruct fi,i = 1, . . . , N at the same time. Figure ^ shows two examples of 
applications in MR images. 

This paper is organized as follows : In section 2 we introduce the common feature 
z, model the relation between the images fi to it through p(fi\z) and its proper 
characteristics through a prior law p(z). In section 3 we give detailed expressions 
of the a posteriori law and propose a general structure of the MCMC algorithm 
to estimate / and z. Finally, in section 4 we present some simulation results to 
show the performances of the proposed method. 
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2 Modeling for Bayesian data fusion 

In this paper we consider the model (2) where after discretization and using the 
notations g t = [g 4 (l), . . .,gi(S)] T and g = (gi) l=1 ....M , ft = [/»(1), • • • , fi{S)] T 
and / = (/i)j=i,...,M, and e, = [e»(l), . . . , ei(S)] T and e = (ei)i=i,...,M, with 5" 
the total number of pixels of the images fi, we have : 

£ = / + £, or gi = fi + £i, i = l,...,M (3) 

Within this model and assuming Gaussian independant noises, p(£i) = A/"(0, erf. J), 
we have 

M M 

p(£i/) = n^i^) = Hpe*tei - fi) 

i=l i=l 

As we want to reconstruct an image with statistically homogeneous regions, it 
is natural to introduce a hidden variable z = (z(l), . . . , z(S)) 6 {1, . . . , K} s 
which represents a common classification of the images fa. The problem is now 
to estimate the set of variables (f,z) using the Bayesian approach : 

P(f,z\g)=p(t\z,g)p(z\g) (4) 

Thus to be able to give an expression for p(f, z\g) using the Bayes formula, we 
need to define p(gi\fi), p{fi\z), p{gi\z) and p(z). 

Assuming £j centered, white and Gaussian, ans S the number of pixels of an 
image, we have : 

p{gi\fi)=M{fi,°li) = (2^5-) 2 cxp|-^|| Sl - /,|| 2 | 

To assign p(fi\z) we first define the sets of pixels which are in the same class : 

Rk = {r : z(r) = k}, \R k \ = n k 
f ik = {Mr) : z(r) = k} 

Then we assume that all the pixels of an image fi which are in the same class 
will be characterized by a mean nn k and a variance of fc : 

P(fi(r)\z{r) = k) =N{m lkl a 2 lk ) 

With these notations we have : 

p{f lk )=M(m lk l,a 2 tk I) 

K 



p(fi\z) = l[Af(m ik l,a* k I) 



fe=i 

K 



II — ^ cxp{--^||/ lfc -m lfc l|| 2 j, i = l,...,M. 



fe =i \ \/z™hJ 1 2<Tik 
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The next step is to define p(gi\z). To do this we may use the relation (3) and 
the laws p(fi\z) and p(£i) to obtain 

p{g l {r)\z{r) = k) = M(m ik , of fc + a 2 e .) 

Finally we have to assign p(z). As we introduced the hidden variable z for finding 
statistically homogeneous regions in images, it is natural to define a spatial 
dependency on these labels. The simplest model to account for this desired local 
spatial dependency is a Potts Markov Random Field model : 

P(z) = exp < a S ( z ^ ~ Z ( S W 

^ ' [ reSseV(r) 

where S is the set of pixels, 5(0) — 1, S(t) = if t ^ 0, V(r) denotes the 
neighborhood of the pixel r (here we consider a neighborhood of 4 pixels) and 
a represents the degree of the spatial dependency of the variable z. This pa- 
rameter will be fixed for our algorithm. We have now all the necessary prior 
laws p(gi\fi), p(fi\z), p{g%\z) and p(z) and then we can give an expression for 
p(f,z\g). However these probability laws have in general unknown parameters 
such as <j\. in p(gi\fi) or m ik and of fe in p(fi\z). In a full Bayesian approach, 
we have to assign prior laws to these "hyperparameters" . Then the choice of 
prior laws for the hyperparameters is still an open problem. In [S] the authors 
used differential geometry tools to construct particular priors which contain as 
particular case the entropic and conjugate priors. In this paper we choose this 
last one. 

Let rrii = (rrii k )k=i,...,K an d = (°f fc)fc=i,...,K be the means and the variances 
of the pixels in different regions of the images fi as defined before. We define 9i 
as the set of all the parameters which must be estimated : 

l = (a\.,mi,(T\), i = l,...,M 

and we note = (6i)i=i,....M ■ When applied the particular priors of ([Hj) for our 
case, we find the following conjugate priors : 



— Inverse Gamma TGia^ > fto' ) anc ^ 2£/(ai0) 0io) respectively for the variances 
a\ % and af k , 

— Gaussian Af(miQ,af ) for the means rrijfc. 

The hyper-hyperparameters ai , f3i Q , mj and of are fixed and the results 
are not in general too sensitive to their exact values. However in case of noisy 
images we can constrain small value on af in order to force the reconstruction 
of homogeneous regions. 

3 A posteriori distributions for the Gibbs algorithm 

The Bayesian approach consists now to estimate the whole set of variables 
(/, z,0) following the joint a posteriori distribution p(f, z,9\g). It is difficult to 
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simulate a joint sample (/, z, 0) directly from his joint a posteriori distribution. 
However we can note that considering the prior laws defined before, we are able 
to simulate the conditionnal a posteriori laws p(f, z\g, 0) and p(9\g, /, z). That 
is why we propose a Gibbs algorithm to estimate (/, z,9), splitting first this set 
of variables into two subsets, (/, z) and (0) : 

p(f,z\g,0)=p(f\z,g,9)p(z\g,0) 

Then the sampling of this joint distribution is obtained by sampling first p(z \g,0) 
and then sampling p(f\z, g, 0_). We will now define the conditionnal a posteriori 
distribution we use for the Gibbs algorithm. 

Sampling z\g,9 : 
for this step we have : 

M 

P(z\g,0) ocp(g\z,0) p(z) = J|p(9i|z,0i) p(z) 

i=i 

As we choosed a Potts Markov Random Field model for the labels, we may note 
that an exact sampling of the a posteriori distribution p(z\g,9) is impossible. 
However we propose in section 4 a parallel implementation of the Gibbs sam- 
pling for resolving this problem. 

Sampling fi\gi,z,0i : 

We can write the a posteriori law p{fi(r)\g i (r) 1 z{r), 0i) as follows : 
p(Mr)\ 9l (r),z(r) = k, 0;) = ^(m^V? T*) 

where 

apost 2 a P° st ( 9i\ r ) . m ik \ , 2 a P ost f 1 . 1 

™i k =^k [—2- + -T J and °ik =[-T + -2 

\ °£, "i fc/ \"si a i k 

sampling 0i\fi,gi,z : 

We have the following relation : 

P(0i\fi,9i,z) ^p{o- 2 e t \ft,gi) p{mi,<r1\fi,z) 

and using again the Bayes formula, the a posteriori distributions are calculated 
from the prior selection fixed before and we have 
- rn ik \fuz,o-1 k ,m i0 ,o-? ~-A/"(/z ife ,t>? fe ), with 

mk = v 2 ik ( 5^ + 4- J2 £( r )) and v ^k=[^r + -j- 

\"i0 ik reRk J \ a ik a i( 

~ °ik\fi> Z ' a iO,0iO ~TG(aik,Pik)> with 

aik = aao + Y and (3 lk = f3 lQ + ^ ^ (/,(r) - m ik f 
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- ofj/i) 9i ~ -Si), with 

5 1 
v% = — + osq*, S 1 = number of pixels and ^ = - /i|| 2 + /3g 4 (5) 

4 Parallel implementation of the Gibbs algorithm 

As we choosed a first order neighborhood system for the labels, we may also 
note that it is possible to implement the Gibbs algorithm in parallel. Indeed, 
we can decompose the whole set of pixels into two subsets forming a chessboard 
(see figure 2). In this case if we fix the black (respectively white) labels, then 
the white (respectively black) labels become independant. This decomposition 




reduces the complexity of the Gibbs algorithm because we can simulate the whole 
set of labels in only two steps. The Parallel Gibbs algorithm we implemented is 
then the following : given an initial state (0%, 02, z)(°\ 

Parallel Gibbs sampling 

repeat until convergence 

1. simulate z B {n) ~p(*|«V (n ~ 1) ,£ s fl (n_1) ) 
simulate zw^ ~ p(z\z"b^ ,g,§. ^) 
simulate ~ p{fi\gi,z^ n \d} n 1} ) 

2. simulate 0^ ~ i (n) ,gi) 

5 Simulation and results 

Here we illustrate two examples of MRI images : PD, Tl-weightd and T2- 
weighted slices of a MR brain image, which are (188 x 193) images for the first ex- 
ample, and Tl-weighted,T2-weightedand Tl-weighted with contrast agent slices 
of a MR brain image, which are (289 x 236) images for the second example. In 
this last we have added a gaussian noise. 

Figures |2| and ^ show the data fusion result of the proposed method. As it is seen 
on these figures the fusionned segmentations we obtain contain all the regions 
and boundaries of the observations, but we have not yet compared with other 
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methods to see the performances of our algorithm. Also the presence of noise in 
figure 0] do not really affect the result of segmentation and, at the same time, 
the proposed algorithm give visibly improved reconstructed images. 




In both applications we have satisfactory results of image fusion, even when 
images present a great number of homogeneous regions and boundaries. Note 
also that in both applications we fixed a prior small value of of to improve 
the reconstructed images. Another way may be the introduction of some local 
spatial dependency between the neighboring pixels of images fi(r). This point 
is under development and we will report soon on the results. 



8 



6 Conclusion 

We proposed a Bayesian method for data fusion of images, with a Potts Markov 
Random Field model on the hidden variable z. We illustrated how a joint seg- 
mentation and reconstruction can be obtained in case of MRI images. We showed 
then how reconstruction and fusion can be computed at the same time using a 
MCMC algorithm. We considered the case of noisy images and showed that the 
joint segmentation is not greatly affected. This method give an unsupervised seg- 
mentation which do not take into account particular shapes and then can give 
good results in many different cases. However we assume for the moment that 
the observed images are registered. We think that this modclisation is promising 
for introducing registration operators Hi and then implementing common seg- 
mentation and registration at the same time. Another perspective is to introduce 
spatial dependency directly on the images /j for involving the reconstruction. 
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